Phase lag in epidemics on a network of cities 
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We study the synchronisation and phase-lag of fluctuations in the number of infected individuals 
in a network of cities between which individuals commute. The frequency and amplitude of these 
oscillations is known to be very well captured by the van Kampen system-size expansion, and we 
use this approximation to compute the complex coherence function that describes their correlation. 
We find that, if the infection rate differs from city to city and the coupling between them is not too 
strong, these oscillations are synchronised with a well defined phase lag between cities. The analytic 
description of the effect is shown to be in good agreement with the results of stochastic simulations 
for realistic population sizes. 
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I. INTRODUCTION 

The theory of the frequency and amplitude of stochas- 
tic oscillations in models of epidemics of childhood dis- 
eases has been extensively developed over the last few 
years [ll-[To|. but the question of the synchrony of these 
oscillations has received comparatively little attention. 
This is despite its undoubted importance; whether the 
oscillations in different locations are in phase or out of 
phase with each other will clearly have consequences for 
the duration of an epidemic and for the persistence of a 
disease. The particular case of in-phase and antiphase 
locking of epidemics in coupled cities has been described 
and interpreted in the literature [ill 113 j but only for 
models that exhibit periodic oscillations in the infinite 
population limit. Similarly, synchronisation in models 
of coupled demographic oscillators has been studied [l^ , 
but not in the context of epidemics. 

This subject has been so little investigated from a the- 
oretical point of view, that it is useful to recall a simpler, 
and more familiar example. It is frequently argued in 
connection with predator-prey interactions that if preda- 
tor numbers happen to be high, prey numbers will sub- 
sequently fall due to increased predation, and when prey 
numbers fall, predator numbers will subsequently fall due 
to lack of food. This is actually an argument which re- 
lates to fiuctuations in the number of predators and prey 
and would seem to indicate oscillations, with predator 
and prey numbers being out of phase with each other. 

A similar argument, although not so well-known, can 
be applied to infected and susceptible individuals in a 
city. In both cases, the frequency and amplitude of these 
oscillations is very well captured by the van Kampen 
system-size expansion [H-l^, fT3 . [T5| , and it seems reason- 
able to suppose that the same approach should also be 
able to quantify the synchrony and phase-lag between os- 
cillations of different kinds of individuals. In this paper 
we apply this method to study possible synchronisation 
between fluctuations in the number of infected individu- 
als in a network of cities between which individuals com- 



mute. 

A previous study [16] has shown that if the parameter 
setting the infection rate, denoted by /?, is the same for 
all cities in the network, then the stochastic dynamics 
takes a remarkably simple form. This has been known 
for some time in the case of the deterministic dynam- 
ics: a rather general theorem tells us that in this case 
the fraction of susceptible, infected and recovered indi- 
viduals is the same in all cities fl?']. In Ref. [T^], it was 
shown that the stochastic oscillations in this situation 
are also quite simple, having only one frequency, which 
is independent of the city, even if the sizes of the cities 
are different. As we will show later, these oscillations 
are synchronised between cities, but with no phase lag, 
in agreement with the findings reported in the literature 
for stochastic simulations of infection dynamics on cou- 
pled population patches [ill, [3 • Exploring spatial het- 
erogeneity as a means to overcome some of the shortcom- 
ings of the standard description of epidemics, namely in 
the predicted frequency of disease fade-outs, was one of 
the goals of these studies. In this sense, the trivial phase 
relation between cities is a negative result. 

However, it was assumed in previous work that the 
infection rate /3 is the same in all cities. If /3 is taken to 
be different in different cities, the symmetry of the fixed 
point is broken, which introduces the possibility of more 
complex oscillations. There are other ways of achieving 
this, but within the framework in which we will work, this 
is a natural way of introducing the symmetry breaking. 
We shall see that, under these conditions, the numbers of 
infected individuals fluctuate with a phase lag between 
cities, provided that the coupling is not too strong. 

The paper is structured as follows. In Sec. |TT] we de- 
scribe the stochastic model, including a microscopic in- 
terpretation of the parameters in terms of the mobility 
patterns of the populations and the characteristics of the 
disease. The deterministic equations in the infinite pop- 
ulation limit are derived in Sec. Illl[ as are the Langevin 
equations for the fluctuations around the non-trivial de- 
terministic equilibrium in the linear noise approximation. 



2 



found using van Kampen's system-size expansion. In 
Sec. lIVI we introduce the complex coherence function that 
measures the cross-correlations between different cities 
and/or types of individuals, and we compute it for sus- 
ceptibility and infection in the simple case of one city as 
an illustration of the method. The correlation between 
infection in two cities is studied in Sec. IVland it is shown 
that for moderate coupling and different /3 the fluctu- 
ations have a characteristic frequency range and phase 
relation. This is extended to the case of three cities in 
Sec. Ell We conclude in Sec. EIIl 



II. MODEL 

The model consists of n cities, labeled j = 1, . . . , n con- 
taining individuals who are infected, susceptible or recov- 
ered from a disease. As is common, for mathematical con- 
venience the total number of individuals in a particular 
city, Nj for city j, is fixed. This means that the number 
of recovered individuals in each city may be expressed 
in terms of the number of infected and the number of 
susceptible individuals: Rj = Nj — Sj — Ij, j ^ 1, . . . ,n. 
This reduces the number of degrees of freedom from three 
to two per city, considerably simplifying the analysis. It 
also means that deaths and births are coupled, with the 
random death of an individual giving rise to the birth of 
a susceptible individual, and so providing a fresh victim 
for the disease. 

We will have in mind a particular type of interchange 
of individuals between cities, but we will see that the 
construction which comes from this is sufficiently general 
to encompass a large variety of situations. A concrete 
example of what we have in mind is a set of residential 
areas and one or more urban areas; there would typically 
be a significant fraction of commuters from the former to 
the latter, and fewer from the latter to the former. The 
model will thus allow for a fraction of the individuals 
from city j to commute to city fc, denoted by fkj, with 
the number of non-commuters (residents) in city j being 
1 — fj = 1 — J2k^ fkj- The model will differ from that 
studied in Ref. [l6| in that we will take into account the 
different nature of the city (for instance, residential or 
urban) by allowing the parameter describing the rate of 
infection, /3, to vary from city to city. 

If we denote the state of the system of n cities at a given 
time by a = {Si, Ii, . . . , Sn, In}, then the recovery of an 
individual in city j will occur at a rate jlj and cause 
a transition to the new state a' = {Si, Ii, . . . , Sj, Ij — 
1, . . . , 5'n, /„}. We will write the transition rate for this 
case as 



T{S,,I, 



(1) 



with the initial state on the right and the final state on 
the left. It should be noted that we have only listed only 
the variables in city j in order to lighten the notation. In 
a similar fashion, the death of an infected or recovered 



individual, and the birth of a susceptible, in city j are 
given by 

T{s, + i,i,-i\s„i,) = 

T{Sj + 1, 1,\S,,I,) - ii{N, - S, - /,), (2) 
respectively. 

The network structure manifests itself when construct- 
ing the transition rates induced by infections. The rate of 
infections involving susceptibles from city j and infectives 
from city k will be proportional to Sjik- There will be 
five types of term. In two of them k — j. These are when 
the infective residents in city j infect susceptible residents 
in city j and when infective commuters from city j infect 
susceptible commuters from city j in city £ {£ ^ j). The 
rates for these are respectively /3j(l — fj^Sjlj/Mj and 
j3if'^^SjIj/Mi, where 13 j is the parameter which charac- 
terises the infection rate in city j and where Mj is the 
number of individuals in city j: 



(l-/,)iV,+ V/,™iV, 



(3) 



The other three terms result when the susceptibles 
from city j are infected by individuals from city k 
where cities j and k are different. Then, infective 
commuters can infect resident susceptibles at a rate 
/3j{l — fj)fjkSjIk/Mj, infective residents can infect com- 
muting susceptibles at a rate Pkfkji^ ~ fk)SjIk/Mk and 
finally infective commuters from city k infect suscepti- 
ble commuters from city j in city i (i ^ j, k) at a rate 
PiftjfikSjlk/Mi. These results allow us to write down 
the transition rate for infection as 



T{S,~\,I, + l\S 



(4) 



where 



Me 



fijk - 



{I - fj) fjkNk , Mkj{l-fk)Nk 



E 



Me 



Mk 

, j,k = l,...,n;j ^k. (5) 



Although we have a particular picture of how individu- 
als move between cites, and of the assignment of infection 
rates to the cities themselves, the final form of the tran- 
sition rate (U) is quite general. Our results will therefore 
apply to a wide range of the possible types of interchanges 
of individuals and choice of infection parameters. Unlike 
the case where the f3j were all equal [1^1 , the (3jk have no 
relations between them, and so in general are indepen- 
dent. Therefore they can be chosen from a consideration 
of Eq. ([SJ or simply externally imposed. 
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III. DYNAMICS 

The model defined in Section is stochastic and 
Markovian, since the transition rates ([T]), ^ and ^ do 
not depend on past states of the system. This means 
that adopting a continuous time description, the time 
evolution of the system may be obtained from a mas- 
ter equation for P{a,t), the probability distribution for 
finding the system in state a at time t [l^ [l^, HO] 



dP{a, t) 
di 



J2 [na\a')P{a',t)^T{a'\a)P{a,t)], (6) 

where the T{aW) are each of the rates ^ and (g]) 
taken in turn [H [11, [IS]. 

The full master equation ^ cannot be solved exactly, 
but the aspects that concern us here can be analysed in a 
remarkably precise way using the system-size expansion 
of van Kampen ■ We have carried out this calculation 
on a simpler version of this model — when all the /3j 
were equal — elsewhere [l^, and while the predictions 
are different, the method of carrying out the expansion is 
very similar. A comparison of the transition rates shows 
that the results for the model we are investigating here 
may be obtained from those in Ref. [l^ by substituting 
Pcjk by Pjk- We will therefore only briefly indicate the 
steps in the analysis, and refer the reader to Ref. [1S| for 
more details. 

The method begins by making the following ansatz 



.J, (7) 

where j = Here Sj = limjv -s-cx) S'j/A'j and 

ij = limAr^._j.oci Ij /Nj are the fraction of individuals from 
city j which are respectively susceptible and infected in 
the deterministic limit. Therefore the variables are bro- 
ken down into a sum of deterministic terms, Sj and ij, 
and the stochastic deviations from these, Xj and j/j . The 
deterministic terms satisfy the ordinary differential equa- 
tions 



dt 
di 



= - V/3jfcSjifc+/i(l-Sj): 



k=l 



^ = '^PjkSjik - {-f + fi)ij, 
fe=i 



(8) 



where j — I, . . . ,n. 

The stochastic fluctuations, Xj and yk, which describe 
the linear fluctuations around trajectories of the de- 
terministic set of equations ([U, are found to obey a 
set of linear stochastic differential equations. For con- 
venience we introduce the vector of these fluctuations 
z — {xi, . . . , Xn, 2/1, ■ • ■ , J/n) and indices J, K = 1, . . . , 2n, 
Then these stochastic differential equations have the 
form di [H iO] 



-^=Y,AjKZK+Vj{t), J=l,...,2n, (9) 



2n 

E 

K=l 



where i]j{t) are Gaussian noise terms with zero mean 
which satisfy {vjit)VKit')) = BjK5{t - t'). The two 
2n X 2n matrices A and B are obtained from the ex- 
pansion procedure. In fact, to determine the matrix A 
it is not necessary to carry out the full system-size ex- 
pansion, since it is related to the Jacobian, , found 
from linear stability analysis about a fixed point of 
Eq. dH]). The precise relationship is J' = S^^AS, where 
S = diag(v^, . . . , v^) [HI. The explicit forms for J 
and B are most easily given in terms of the four n x n 
submatrices: 



J 



Jil) 




Ji3) 


^(4)_ 



(10) 



and similarly for B. The elements of these submatrices 
are found to be 



d jk ^ ''Jh'jki 
n 



Jjk = -Ui + l)5jk + SjPjk, 



and 



(11) 



r(2) _ r,(3) 

^3k -Bjk 



-fiSjkij - 5jk Sj 



Bjk = {7 + fi)Sjkij+SjkYsjl3jiie. (12) 

These matrices depend on the solutions Sj and ij of 
Eq. ([8]), and so on time, but since we will be interested in 
fluctuations about the stationary state, we use the fixed 
point values of Sj and ij, denoted as s* and i*. The re- 
sulting matrices, J"* and B* , are then time- independent. 



IV. ANALYSIS 

Adding the two sets of equations in ([8]), we immedi- 
ately see that the fixed points satisfy 

{j + fi)t*=fi{l-s*), j = l,...,n. (13) 

Using this equation to eliminate the i*, one finds that 



(7 + M) + E^:'fc(l~^* 



fe=i 



= (7 + /-i), j = l, 



(14) 

We will assume that the matrix of the coupling co- 
efficients f3jk is irreducible, which means that any two 
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cities have a direct or indirect interaction. Otherwise the 
n cities may be spUt into non-interacting subsets, and 
several equiUbria may be found by combining disease ex- 
tinction in some subsets with non-trivial equilibrium in 
other subsets. The disease extinction fixed point is sim- 
ple to find. If we assume that any one of the j* is zero, 
for instance = 0, then from Eq. s| = 1. From 
Eq. dS]) we see immediately that J2k=i f^e^ik = 0- Since 
the matrix of couplings is irreducible, and from Eq. ^ 
the entries are non-negative, it follows that — for all 
k. 

Although it is difficult, and in most cases impossible, 
to determine any non-trivial fixed points analytically, a 
theorem tells us that a unique stable fixed point ex- 
ists and is stable. The linear stochastic deviations from 
this fixed point satisfy Eq. © , with A and B being evalu- 
ated at the fixed point. To find the dominant frequencies 
of the stochastic oscillations, and as we shall see also to 
study synchronisation and phase-lag, we Fourier trans- 
form Eq. ^ to obtain 



{-iujSjK - Ajk) zk{lo) = ryj(a;), J = 1, . . . , 2n, 

K=l 

(15) 

where the / denotes the Fourier transform of the function 
/. Defining the matrix —iujSjk — Ajk to be ^jif (w), the 
solution to Eq. is 



2ti 
K=l 

We now introduce the matrix 



(16) 



Pjk{(^) = {z.j{uj)z*k{uj)) 

2n 2n 

= 1111^:jI^^)Blm{^%I^H. (17) 

L=l M=l 

where * now denotes complex conjugation. 

In previous studies, where the focus was on finding 
the frequencies and amplitudes of the stochastic oscilla- 
tions jlnSi, tl^] , only the power spectrum (when J = K) 
was analysed. Here we will also be interested in the cross- 
correlations between infection in two different cities, and 
so will also wish to calculate the cross-spectrum (when 
J ^ A'). It is frequently convenient to normalise this by 
the relevant power-spectrum, and instead work with the 
complex coherence function (CCF) defined by [2l] . [2^ 



C.,K{i0) 



PjKiio) 



^Pjj{uj)PKK{i^) 



(18) 



The CCF will in general be complex for J ^ K , and 
so typically one calculates its magnitude and phase, that 
is, the coherence 



\Cjk{u^)\ 



\Pjk{uj) 



y/Pjj{uj)PKK{uj)' 



(19) 
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FIG. 1: (Color online) Synchronisation and phase-lag between 
susceptible and infected individuals in the one-city SIR model. 
Top: Analytical power spectra of the fluctuations of suscepti- 
bles (black line), Pii(a;), and of infectives [cyan (gray) line], 
^"22(1^), given by Eq. (|17l) . The dotted vertical line indicates 
the frequency Umax at which these spectra attain maxima. 
Bottom: Phase spectrum (left), (^21(1^), and coherence (right), 
|C2i(a;)|, for the fluctuations of infectives and susceptibles. 
The solid lines are the analytical results given by Eq. (|22p 
and Eq. (|2ip . The open circles are the results for the same 
quantities obtained from simulations. Parameters: TV = 10®, 
/? = 17(7 + fj,), n = 1/50 1/y and 7 = 365/13 1/y In all 
panels the dashed vertical lines bound the frequency range 
where the stochastic amplification is significant. 



and the phase spectrum 

'ImiCjKiuj)) 



(t).JK{<^) = tan" 



Re(CjK(w)) 



tan 



Im {P.jk{^)) 



Re(PjK(w))J 
(20) 

As an example of using the CCF to understand syn- 
chronisation and phase-lag in systems with sustained 
stochastic cycles, we apply it to the SIR model in a single 
city. We could also apply it to the predator-prey system 
discussed earlier, but we already have the required equa- 
tions for the one-city SIR model in this paper: Eqs. ([8])- 
(ITil) . with the indices omitted and with /J^-j. replaced by /3. 
It should be emphasised that in this example we are look- 
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ing at the synchronisation and phase-lag between suscep- 
tible and infected individuals, whereas in the actual ap- 
plication of the formalism to n cities, the interest is in 
possible synchronisation and phase-lag between infected 
individuals in city j and infected individuals in city k. 

For models such as the one-city SIR and the predator- 
prey model, Eq. ([T9| becomes 



IC21MI 



' g4Uj'^ + ff2^^ + go 

h4UJ^ + h2Uj^ + ho ' 



(21) 



where the coefficients of the polynomials are sums of 
products of the Ajk and Bjk, J,K — 1,2, which can be 
found from Eqs. P^ - (IT^ . and Eq. ([20)1 becomes simply 



021(1^) = tan 



ko + k2U}'^ 



(22) 



where again ko , fci and /c2 are sums of products of the 
Ajk and Bjk, J,K = 1,2. 

The power spectra of the fluctuations of susceptible 
and infected individuals computed from Eqs. (fTTl) . Pii(w) 
and P22(w), are shown in the top panel of Fig. 1. Stochas- 
tic amplification gives rise to significant fluctuations in a 
well-defined frequency range centered at the frequency 
'-^max where both Pii{uj) and P22(w) peak. The coher- 
ence, |C2i(a;)|, and phase spectrum, (j)2i{uj), given by 
Eq. (PT)) and Eq. (1^ are shown in the right and left 
bottom panels of Fig. 1, together with the results for 
the same quantities obtained from numerical simulations. 
There is very good agreement between the analytic ap- 
proximation and the simulations in the frequency range 
where the fluctuations have significant amplitude. Out- 
side of this range, the magnitudes of the power spec- 
tra are so low that errors become significant, leading to 
agreement which is not as good. However within this 
frequency range the fluctuations of the two kinds of in- 
dividuals are strongly correlated and have a well-defined 
phase lag. 

The numerical results presented here and in the fol- 
lowing sections were obtained from long numerical sim- 
ulations based on the Gillespie algorithm [2^. Each run 
started from a random initial condition, with the vec- 
tor of the fluctuations of susceptibles and infectives in n 
cities, z — (xi, . . . , Xn, yi, . . . , j/„), computed from sim- 
ulated time series according to Eq. ([7]), and recorded 
at equal time intervals. From each simulation run the 
power and cross spectra for the fluctuations are com- 
puted as zj{uj)z'^{uj), where J, K = 1, . . . , 2n and where 
* denotes complex conjugation, by using discrete Fourier 
transforms. The final spectra are averages of 10'^ to 
5 X 10'^ simulations. Having computed these, the nu- 
merical CCFs, coherence and phase spectra can easily be 
computed from Eqs. p^ - (|20|) . The values for the epi- 
demiological and demographic parameters used in this 
paper are those relevant for measles. In all simulations, 
we take fi = 1/50 1/y and 7 = 365/13 1/y [H-i^. 
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FIG. 2: (Color online) Synchronisation and phase-lag between 
infected individuals in city 1 and infected individuals in city 2 
in the two-city SIR model. Top: Power spectra of the fluctua- 
tions of infectives in city 1, P33(ij), and of infectives in city 2, 
P4,4,{td), given by Eq. (|17[). The dotted vertical lines indicate 
the frequencies Uma^i and u}max2 at which these spectra at- 
tain maxima. The dashed vertical lines bound the frequency 
range where the stochastic amplification is significant. Bot- 
tom: Parametric plot of the CCF, (743(0;), in the complex 
plane. The solid line is the analytical result given by Eq. psp 
[the parts of the curve shown in red (black) and in gray cor- 
respond to 1 < a; < 5 and {0 < < 1} U {5 < w < 100}, 
respectively] and the open circles are the result for the same 
quantity obtained from simulations for < w < 10. Param- 
eters: ^1 = 12(7 + fi), 132 = 17(7 + ^i), A^i = iVa = 10^ 
/12 = 0.001, /21 = 0.01, fj. = 1/50 1/y and 7 = 365/13 1/y 



V. TWO CITIES 

We now turn to the description of the synchronisation 
and phase-lag between infected individuals in two dif- 
ferent cities, using the quantities introduced in Sec. IIVI 
Plots of the power spectra for infectives, P33(a;) and 
P44(w), are shown in the top panel of Fig. 2 for a cer- 
tain choice of parameters (recall that Z3 and are the 
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FIG. 3: Dependence of the synchronisation and phase-lag 
between infected individuals in city 1 and infected individ- 
uals in city 2 on the infection rates in the two-city SIR 
model. Top: Phase (left), (jfisi^M), and coherence (right), 
|C43(i^Af)|, given by Eq. (f20)) and by Eq. (fT9)l as a function of 
(/32 — Pi) / {"f + jj.) . Bottom: The difference between the peak 
frequencies ujmaxi of P33 and iJmax2 of Pn as a function of 
(/32 — /3i ) / (7 + /i) • To obtain the plots /32 / (7 + /^) was fixed at 
15 and I3i/{'y + fi) was varied between 5 and 32.5. Parameters: 
iVi = 10", N2 = 2x 10^ fi2 = 0.005, /21 = 0.01, fi = 1/50 
1/y and 7 = 365/13 1/y. 



fluctuations of infection in cities 1 and 2 respectively). 
Different values have been taken for /3i and /32 , to reflect 
different social contact patterns in the two cities. In con- 
trast with the spectra for different types of individuals 
in the one-city case, P33(ci;) and P44(a;) peak at different 
frequencies ujmaxi and U!max2 ■ The frequency range of in- 
terest is bounded by the two dashed vertical lines in the 
figure — outside this range, P33 and P44 have negligible 
amplitude. A parametric plot of the range of the coher- 
ence function (743(0;) in the complex plane is shown in 
the bottom panel, where the darker (red in the color ver- 
sion) portion of the full line corresponds to the frequency 
range highlighted in the top panel. It can be seen that the 
fluctuations in the number of those infected in the two 
cities are correlated and synchronize with a well-defined 
phase relation in the whole frequency range where both 
their amplitudes are significant. 

Also shown in the bottom panel of Fig. 2 are the results 
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FIG. 4: Dependence of the synchronisation and phase-lag be- 
tween infected individuals in city 1 and infected individuals 
in city 2 on the fractions of commuters between the cities 
in the two-city SIR model. The description and coloring are 
as in Fig. 3. To obtain these log-linear plots, /12 was fixed 
at 0.001 and /21 was varied between and 1. Parameters: 

i3i = 12(7 + fi), /?2 = 17(7 + 11), Ni = lo^ N2 = 2x lo^ 

/i = 1/50 1/y and 7 = 365/13 1/y. 



for Ci3{uj) obtained from numerical simulations. There 
is a good agreement between simulations and analytic re- 
sults within the frequency resolution limits of the former, 
showing that for the chosen system sizes the fluctuation 
cross-correlation behavior is well captured by the linear 
noise approximation. 

In order to investigate the dependence of the coher- 
ence and phase spectrum on the choice of parameters, 
we have computed the value ujm for which |C43(a;)| at- 
tains its maximum, |C43(cjAf)| and 043 (wm)- A plot of 
(/>43(wA/) and |C43(i:^Af)| as a function of (/32 - /3i)/(7 + A^) 
is shown in Fig. 3, top panels. To obtain the plot only 
the infection rate (32 was varied while the remaining pa- 
rameters were kept fixed. With respect to the example of 
Fig. 2, we have chosen one of the populations, N2, to be 
twice as large and one of the coupling parameters, /12, 
to be five times larger. It should be noted that when 
{(32 — /3i ) / (7 + ^) is so large that the frequency ranges 
of the stochastic fluctuation peaks in each city no longer 
overlap, |C43(w)| becomes bimodal instead of unimodal 
and ujm is no longer well defined. The difference be- 
tween the peak frequencies ujmaxi of P33 and uJmax2 of 
P44 is shown in the bottom panel of Fig. 3. For values 
of {(32 — /?i) /(7 + ^) to the left of the represented range, 
i.e. for values smaller than —17.5, this difference becomes 
larger leading to a more complex synchronization pattern 
and a bimodal |C43(w)|. 

We also see that, in this case, fluctuations with a non- 
trivial phase relation require the two cities to have dif- 
ferent values of the infection rates (3i and (32- We have 
found this to be true in general. More precisely, the imag- 
inary parts of the cross spectra of infection Pj/f (w) with 
J,K — n+1, . . . ,2n and J ^ K, as given by Eq. (fTTl) . are 
identically zero when the infection rates « = 1, . . . , n. 
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are the same in all cities, independently of the choice 
of the population sizes and fractions of commuters. The 
proof of this result is given in the Appendix. This is com- 
patible with either in-phase or antiphase synchronization. 

Another condition for (^43 (wm) to be different from 
zero is that the coupling between the two cities is not too 
strong. A plot of |C43(wAf)| and (^43(1^ a/) as a function 
of the coupling /21 is shown in Fig. 4. The remaining pa- 
rameters are as in Fig. 2, except for N2, which was taken 
twice as large. It can be seen that |C43(cjj\f)| increases 
monotonically with /21, indicating an increase in cross 
correlation between infected individuals in the two cities 
as the coupling gets stronger, but 043(0; a/) tends to zero 
with increasing /21. 



VI. THREE CITIES 

The analysis carried out in Sec. |V] can be extended to 
three or more cities. In this Section we will illustrate 
this by investigating two examples for three cities. We 
have checked that the approximate analytic description is 
again in good agreement with the results of simulations, 
and the plots show the analytic results only, to avoid 
cluttering. 

In the case of Fig. 5, the population sizes are equal, the 
demographic coupling parameters fij are all equal, but 
the three transmission rates j3j, j = 1,2,3 are different. 
The power spectra of the fluctuations of infectives in each 
city, Pjj, J = 4, 5, 6, shown in the top panel of the figure, 
behave as in the example of Fig. 2, as do the coherence 
and the phase spectra. The ranges of C54, Ce4 and Cgs 
in the complex plane are plotted in the bottom panel of 
Fig. 5. The phase lags between cities 4 and 5 and 5 and 
6 are approximately the same in this case. 

In the example of Fig. 6 we have considered differ- 
ent population sizes and coupling parameters /y and we 
have also increased the difference between the transmis- 
sion rates /3j, j — 1,2,3 with respect to the previous 
example. The results for the fluctuation power spectra 
and for the ranges of C54, Cq^ and Ces in the complex 
plane are plotted in top and bottom panels of the Fig- 
ure. In this case we find three different phase lags for the 
correlated fluctuations in the three city pairs. 

These phase lag effects have implications for estimates 
of the duration of an epidemic and of the likelihood of dis- 
ease extinction. In a metapopulation model with differ- 
ent transmission rates epidemic bursts should last longer, 
and disease extinction should occur less frequently than 
in a single population with the same overall size. 



VII. CONCLUSIONS 

In this paper we have considered a stochastic metapop- 
ulation version of a susceptible-infected-recovered model 
representing infection dynamics in different demograph- 
ically coupled urban centers. We derived the state tran- 
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FIG. 5: Synchronisation and phase-lag between infected in- 
dividuals in three cities with different infection rates in the 
three-city SIR model. Top: Power spectra of the fluctuations 
of infectives in city 1, P44{uj), in city 2, Ps5{iu), and in city 3, 
P66(i^), given by Eq. (|17p . The dotted and the dashed vertical 
lines have the same meaning as in Fig. 2. Bottom: Paramet- 
ric plot of the CCFs given by Eq. 1181 C54(tj), Cei{Ld), and 
Ce,5{ijj), in the complex plane. Parameters: j3i = 7(7 + n), 
Pi - P2. P3 = 7: 12: 17, iV> = 10^ ^ = 1/50 1/y, 7 = 365/13 
1/y, and fij = 0.005, where i = 1, 2, 3 and i / j. 



sition rates of the stochastic process from a microscopic 
model for the mobility of individuals between cities. Sim- 
ilar rate equations based on a phenomenological coupling 
parameter have been proposed in the literature. Here we 
have adopted a bottom-up approach, as an illustration of 
how the parameters of the general model can be related 
with the mobility patterns of the sub-populations. 

The correlations between the fluctuations in the num- 
ber of infected in different cities were studied by means 
of the complex coherence function. It was found that, if 
the infection rate differed from city to city and the cou- 
pling was not too strong, oscillations were observed which 
were synchronised with a well defined phase lag between 
cities. We also showed that for realistic population sizes 
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CO 



on daily monitoring of urban populations . The con- 
ditions found in this paper for a non-trivial phase re- 
lation between the disease incidence fluctuations in con- 
nected urban centres can therefore be considered realistic 
in many settings. 
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FIG. 6: An example of synchronisation and phase- lag between 
infected individuals in three cities in the case when all param- 
eter values are different. The description and coloring are as 
in Fig. 5. Parameters: /3i = 7(7 -h/x), /3i : /32 : /Js = 7 : 12 ; 27, 
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7 = 365/13 1/y, /12 = /31 = 0.001, /13 = hi = 0.002, 
/21 = 0.01, and /32 = 0.005. 



this effect was well described analytically by the linear 
noise approximation. 

The combined effect of stochasticity and demographic 
coupling as a possible driver of correlations and spatial 
patterns that are missed by traditional epidemic models 
was suggested long ago O, [13 • This idea has been left 
largely unexplored because these early studies based on 
computer simulations took the infection rate to be the 
same in different cities. In this case, as we have seen, the 
fluctuations are trivially synchronized, with no phase lag. 
However, the infection rate is a phenomenological param- 
eter that depends not only on the disease but also on the 
rate of potentially infectious contacts that characterize 
a given social environment. In particular, a strong pos- 
itive correlation between measles transmission rate and 
population density was confirmed in recent study based 



Appendix A: Proof that P*^^' is real when the 
infection rates are equal 

In this Appendix we will show that the matrix Pjk 
defined by Eq. (ITTl) is real for J,K = n + 1, . . . ,2n in 
the case where the infection rates for each city are equal: 
Pj — /3 for all j. This shows that the cross-spectra for 
fluctuations in the number of infected individuals is real, 
and so that in this case there is either no phase lag or a 
phase lag of tt. We would expect that the former situa- 
tion, i.e. no phase lag, would hold, and we have explicitly 
checked this for the case of two and three cities. For gen- 
eral n we have a partial proof that the phase lag is zero, 
and work is underway to complete it. 

The results for the situation where /Sj — /3 for all j can 
be inferred from the results of the current paper, and are 
given explicitly in Ref. [l6j . They are given in terms of 
the matrices A and B at the fixed point and can be most 
easily expressed in terms of four n x n submatrices by 
writing 



A 



Ad) 


A(2) 


A(3) 


^(4) 



(Al) 



and similarly for B and = —iuil — A. Here / is the 
n X n identity matrix. The matrices and ^^^^ are 
proportional to the identity matrix and we write them as 
pi and a I respectively. The inverse of then takes 

the simple form 



" _$(4) 


$(2) " 


al 


-Pl 



Y 








Y 



(A2) 



where is the n x n matrix (7$ (2) _ 

In the calculation of Pjxi'^) the matrices Y appear in 
the combination 



Z = 



Y 







bW 


i?(2) 





Y 




B(3) 


B(4) 














(A3) 
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However in the case we are interested in [16|, the sub- 
matrices i?*^"), a = 1, ... ,4 are diagonal: -B]"'* = ^a^jk 
with j, fc = 1, . . . , n. Therefore 



(A4) 



BiYY^ 


B2YY^ 




B^YY^ 



by 



The quantity we wish to study, Pjk{(^), is now given 



" _$(4) 


$(2) " 


cr/ 


-pi 



"^(1) 


Z(2)" 


^(3) 





' _$(4)t 


al 


$(2)t 


-pi 



p = 



(A5) 

where Z" is one of the four submatrices given in Eq. (|A4p . 
It follows from Eq, (|A5|) that 



p(4) ^ ^2^(1) _ ^ (p + + |p|2z(4), 



(A6) 



where we have used Z'^'^^ = Z^^^ which follows from the 
the symmetry of the matrix B (B2 = -63), and also used 
the reality of a. We now show that the matrix Z is real, 
and so that P*^"*^) is real. 

To do this we write down the explicit forms for the 
entries of (f>(a;). They are 



M/3 
M + 7 



ILO, 



/i 



M/3 



M + 7 

= (/i + 7-?;a;)/- (^ + 7)X, 



(A7) 



1 /2 

where X is an n X n matrix given by XjTc = (N^/Nk) Cjk, 
and is real and symmetric (using Eq. (B5) of |16|). There- 
fore from the definition of the matrix Y 



with 



Y-' = G{Lu)I + Hiu:)x, 



M + 7 



(A8) 



iuj] + 7 - iuj) (A9) 



and 

H{uj) = (Ai + 7) (Ai - iuj) ■ (AlO) 
Using the symmetry and reality of the matrix 

y-itr-i ^ |GH|2/-H|i/(a,)|V 

+ {G{lo)H{u:) + G{u:)H{lo)]x, (All) 

which is clearly real. Therefore the inverse of this matrix, 
is also real. From Eq. (jA4p it follows that Z is real, 
since B is real. 

We end by showing that the cross-spectra for the fluc- 
tuations in the number of susceptible individuals is also 
real in this case. Although we have not explicitly investi- 
gated this quantity in the main text, it is also of interest, 
and by analogy with P*-^' we would expect it to be real 
too. From Eq, (|X5t 



P(l) = $(4)^(l)$(4)t _ $(4)^(2)^(2)t 

_ $(2)^(3)$(4)t +$(2)^(4)^(2)t^ (A^2) 

Since from Eq. dSl), = {p. + -i -iuj)I -^^^^ and $(2) 
is real, and since the are real, the imaginary part of 
p(i) is given by 

ImP(i) - c^{z(i)co(2)_a>(2)z« 

+ Z(2)$(2)_$(2)Z(3)|. (A13) 

Writing = B^YY"^ and $(2) ^ (^ + -^)^ yields 

ImP(i) = uJ{^i + 7) (Si + P2) { {YY^) X-X (YY^) } , 

(A14) 

where we have used B^ — B2- Now from Eq. (jAll[) . the 
matrix y^^ty-i ig ^ sum of matrices proportional to 
the identity, x s-iid. x^- Therefore it commutes with the 
matrix x, that is, xY~^'^Y~^ ~ Y^^'^Y~^Xi from which 
it follows that YY'^x = xYY\ and so from Eq. (|A14p . 
ImP(i) = 0. 
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